######################
#era5_summary
#paul stainier
#started 2/14/2020
#make summary maps of weather data
#uses old version (2.8) of GADM map https://gadm.org/download_country_v2.html
#map help : https://edzer.github.io/sp/
#####################
rm(list = ls())
library(plyr)
library(dplyr)
library(lubridate)
library(sp)
library(ggplot2)
library(sqldf)
library(RColorBrewer)
library(ggmap)
library(raster)
library(geosphere)
library(forestmangr)
library(tidyverse)
library(data.table)
library(DescTools)
library(hrbrthemes)
library(wesanderson)
library(usmap)
library(ggplot2)
library(mapproj)
library(ggthemes)
library(ggsci)
library(classInt)

#########################
#toggle T/F for the part that 
#you want to run 
#########################
weather_maps <- F
weather_maps_gs <- T

#########################
#replace with own directories before running
#########################
weather_input_directory <- 
analysis_input_directory <- 
output_directory <- 
crosswalk_directory <- 



#########################
#make maps of weather
#for the full year
#########################
if(weather_maps){
  for(weather_var in c("tmax", "prec")){
    setwd(weather_input_directory)
    weather_data <- read.csv(paste("era5", weather_var, "2002_2012_distid3_month_key.csv", sep = "_"))
    if(weather_var == "prec"){
      weather_data$prec = weather_data$prec * 1000
    }
    weather_data <- weather_data[which(weather_data$year <= 2011 & weather_data$year >= 2002),]
    setwd(crosswalk_directory)
    map_distid3_crosswalk <- read.csv("mapdistrict_distid3_crosswalk.csv")
    district_summary <- weather_data %>%
      group_by(distid3, year) %>%
      summarise_all(sum, na.rm = T)
    if(weather_var == "tmax"){
      district_summary$tmax_value <- district_summary$tmax_value / 12
    }
    district_summary <- merge(district_summary, map_distid3_crosswalk, by = "distid3")
    district_summary <- district_summary %>% 
      group_by(district_map_id) %>%
      summarise_all(mean, na.rm = T)
    if(weather_var == "tmax"){
      district_summary$tmax_gt90 <- district_summary$tmax_90_100 + district_summary$tmax_100_110 + district_summary$tmax_gt110
      district_summary$tmax_gt100 <- district_summary$tmax_100_110 + district_summary$tmax_gt110
    }
    india_shapes <- readRDS("IND_adm2.rds")
    #make sure to use ID_2 and not OBJECTID: there is a repeated ID_2 (364 - Palghar)
    #so the numbers for the OBJECTID are off after 364
    colnames(district_summary)[1] <- "ID_2"
    india_shapes <- merge(india_shapes, district_summary, by = "ID_2")
    #tmax value breaks
    if(weather_var == "tmax"){
      breaks.qt <- classIntervals(india_shapes$tmax_value, n = 6, style = "quantile", intervalClosure = "right")
    }
    setwd(paste(output_directory, "maps", sep = "/"))
    if(weather_var == "tmax"){
      pulse_palette <- brewer.pal(n = 7, name = "YlOrRd")
      png("tmax_gt90_map_fullyear_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt90", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_80_90_map_fullyear_2002_2011.png")
      plot(spplot(india_shapes, "tmax_80_90", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_gt100_map_fullyear_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt100", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_100_110_map_fullyear_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_100_110", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_gt110_map_fullyear_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt110", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_value_map_fullyear_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_value", col.regions = pulse_palette, at = breaks.qt$brks, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
    }else if(weather_var == "prec"){
      pulse_palette <- brewer.pal(n = 7, name = "Blues")
      png("yearly_prec_map.png")
      plot(spplot(india_shapes, "prec", col.regions = pulse_palette, cuts = 6,  main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
    }
  }
}

#########################
#make maps of weather
#for the growing season
#########################
if(weather_maps_gs){
  for(weather_var in c("tmax", "prec")){
    setwd(weather_input_directory)
    weather_data <- read.csv(paste("era5", weather_var, "2002_2012_distid3_month_key.csv", sep = "_"))
    if(weather_var == "prec"){
      weather_data$prec = weather_data$prec * 1000
    }
    weather_data <- weather_data[which(weather_data$year <= 2011 & weather_data$year >= 2002 & weather_data$month >= 6),]
    setwd(crosswalk_directory)
    map_distid3_crosswalk <- read.csv("mapdistrict_distid3_crosswalk.csv")
    district_summary <- weather_data %>%
      group_by(distid3, year) %>%
      summarise_all(sum, na.rm = T)
    if(weather_var == "tmax"){
      district_summary$tmax_value <- district_summary$tmax_value / 7
    }
    district_summary <- merge(district_summary, map_distid3_crosswalk, by = "distid3")
    district_summary <- district_summary %>% 
      group_by(district_map_id) %>%
      summarise_all(mean, na.rm = T)
    if(weather_var == "tmax"){
      district_summary$tmax_gt90 <- district_summary$tmax_90_100 + district_summary$tmax_100_110 + district_summary$tmax_gt110
      district_summary$tmax_gt100 <- district_summary$tmax_100_110 + district_summary$tmax_gt110
    }
    india_shapes <- readRDS("IND_adm2.rds")
    #make sure to use ID_2 and not ID_2: there is a repeated ID_2 (364 - Palghar)
    #so the numbers for the ID_2 are off after 364
    colnames(district_summary)[1] <- "ID_2"
    district_summary$tmax_gt100_name = "0"
    district_summary$tmax_gt100_name[which(district_summary$tmax_gt100 > 0 & district_summary$tmax_gt100 <= 10)] = "0-10"
    district_summary$tmax_gt100_name[which(district_summary$tmax_gt100 > 10 & district_summary$tmax_gt100 <= 20)] = "10-20"
    district_summary$tmax_gt100_name[which(district_summary$tmax_gt100 > 20 & district_summary$tmax_gt100 <= 30)] = "20-30"
    district_summary$tmax_gt100_name[which(district_summary$tmax_gt100 > 30 & district_summary$tmax_gt100 <= 40)] = "30-40"
    district_summary$tmax_gt100_name[which(district_summary$tmax_gt100 > 40 & district_summary$tmax_gt100 <= 50)] = "40-50"
    district_summary$tmax_gt100_name[which(district_summary$tmax_gt100 > 50)] = ">50"
    india_shapes <- merge(india_shapes, district_summary, by = "ID_2")
    #tmax value breaks
    if(weather_var == "tmax"){
      breaks.qt <- classIntervals(india_shapes$tmax_value, n = 6, style = "quantile", intervalClosure = "right")
      breaks.qt_gt100 <- classIntervals(india_shapes$tmax_gt100, n = 7, style = "quantile", intervalClosure = "right")
      breaks.qt_gt100$brks <- c(0, 5, 10, 15, 20, 25, 30, 78)
    }

    
    setwd(paste(output_directory, "maps", sep = "/"))
    if(weather_var == "tmax"){
      pulse_palette <- brewer.pal(n = 8, name = "YlOrRd")
      png("tmax_gt90_map_grow_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt90", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_80_90_map_grow_2002_2011.png")
      plot(spplot(india_shapes, "tmax_80_90", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_gt100_map_grow_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt100", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_100_110_map_grow_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_100_110", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_gt110_map_grow_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt110", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      
      png("tmax_gt100names_map_grow_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_gt100_name", col.regions = pulse_palette, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
      

      png("tmax_value_map_grow_2002_2011.png", width = 800, height = 600)
      plot(spplot(india_shapes, "tmax_value", col.regions = pulse_palette, at = breaks.qt$brks, cuts = 6, main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
    }else if(weather_var == "prec"){
      pulse_palette <- brewer.pal(n = 7, name = "Blues")
      png("yearly_prec_map_grow_2002_2011.png")
      plot(spplot(india_shapes, "prec", col.regions = pulse_palette, cuts = 6,  main = "", colorkey = list(labels = list(cex = 3))))
      dev.off()
    }
  }
}

